# Replication for:
# Jesus was a Refugee: Religious Values Framing Can Increase Support for Refugees Among White Evangelical Republicans


# ---------------------------------------------------------------------------------------------------------------
# Appendix Figure 1. Average perceived effectiveness and support of pilot frames
# ---------------------------------------------------------------------------------------------------------------
tgc <- summarySE(Cloud_MTurk, measurevar="human_os", conf.interval = .90,na.rm = T)
tgc$term <- "human"
tgc1 <- summarySE(Cloud_MTurk, measurevar="tol_os", conf.interval = .90,na.rm = T)
tgc1$term <- "tolerance"
tgc2 <- summarySE(Cloud_MTurk, measurevar="moral_os", conf.interval = .90,na.rm = T)
tgc2$term <- "moral"
tgc4 <- summarySE(Cloud_MTurk, measurevar="govt_os", conf.interval = .90,na.rm = T)
tgc4$term <- "gov't"
tgc5 <- summarySE(Cloud_MTurk, measurevar="church_os", conf.interval = .90,na.rm = T)
tgc5$term <- "church"
tgc6 <- summarySE(Cloud_MTurk, measurevar="christian_os", conf.interval = .90,na.rm = T)
tgc6$term <- "christian"
tgc7 <- summarySE(Cloud_MTurk, measurevar="vulnerable_os", conf.interval = .90,na.rm = T)
tgc7$term <- "vulnerable"
tgc8 <- summarySE(Cloud_MTurk, measurevar="policy_os", conf.interval = .90,na.rm = T)
tgc8$term <- "policy"
tgc9 <- summarySE(Cloud_MTurk, measurevar="econrisk_os", conf.interval = .90,na.rm = T)
tgc9$term <- "economic risk"
tgc10 <- summarySE(Cloud_MTurk, measurevar="capacity_os", conf.interval = .90,na.rm = T)
tgc10$term <- "capacity"
tgc11 <- summarySE(Cloud_MTurk, measurevar="muslim_os", conf.interval = .90,na.rm = T)
tgc11$term <- "muslim"
tgc12 <- summarySE(Cloud_MTurk, measurevar="terrorist1_os", conf.interval = .90,na.rm = T)
tgc12$term <- "terrorist 1"
tgc13 <- summarySE(Cloud_MTurk, measurevar="terrorist2_os", conf.interval = .90,na.rm = T)
tgc13$term <- "terrorist 2"
tgc14 <- summarySE(Cloud_MTurk, measurevar="cultural_os", conf.interval = .90,na.rm = T)
tgc14$term <- "cultural"
tgc15 <- summarySE(Cloud_MTurk, measurevar="aid_os", conf.interval = .90,na.rm = T)
tgc15$term <- "aid"
tgc16 <- summarySE(Cloud_MTurk, measurevar="borders_os", conf.interval = .90,na.rm = T)
tgc16$term <- "borders"
tgc17 <- summarySE(Cloud_MTurk, measurevar="disease_os", conf.interval = .90,na.rm = T)
tgc17$term <- "disease"
tgc18 <- summarySE(Cloud_MTurk, measurevar="domore_os", conf.interval = .90,na.rm = T)
tgc18$term <- "do more"
tgc19 <- summarySE(Cloud_MTurk, measurevar="christianreason_os", conf.interval = .90,na.rm = T)
tgc19$term <- "christian reason"
tgc20 <- summarySE(Cloud_MTurk, measurevar="econ_os", conf.interval = .90,na.rm = T)
tgc20$term <- "econ"


Cloud_MTurks <- c("tgc","tgc1","tgc2","tgc4","tgc5","tgc6","tgc7","tgc8","tgc9","tgc10","tgc11","tgc12","tgc13","tgc14","tgc15","tgc16","tgc17","tgc18","tgc19","tgc20")

for(tgcs in Cloud_MTurks)
  assign(tgcs, setNames(get(tgcs),  c("ID", "N", "mean","sd","se","ci","term")))

MASTER <- rbind(tgc,tgc1,tgc2,tgc4,tgc5,tgc6,tgc7,tgc8,tgc9,tgc10,tgc11,tgc12,tgc13,tgc14,tgc15,tgc16,tgc17,tgc18,tgc19,tgc20)

MASTER$low = MASTER$mean-MASTER$ci
MASTER$high = MASTER$mean+MASTER$ci
rm(list=Cloud_MTurks)

tgc <- summarySE(Cloud_MTurk, measurevar="human_eff", conf.interval = .90,na.rm = T)
tgc$term <- "human"
tgc1 <- summarySE(Cloud_MTurk, measurevar="tol_eff", conf.interval = .90,na.rm = T)
tgc1$term <- "tolerance"
tgc2 <- summarySE(Cloud_MTurk, measurevar="moral_eff", conf.interval = .90,na.rm = T)
tgc2$term <- "moral"
tgc4 <- summarySE(Cloud_MTurk, measurevar="govt_eff", conf.interval = .90,na.rm = T)
tgc4$term <- "gov't"
tgc5 <- summarySE(Cloud_MTurk, measurevar="church_eff", conf.interval = .90,na.rm = T)
tgc5$term <- "church"
tgc6 <- summarySE(Cloud_MTurk, measurevar="christian_eff", conf.interval = .90,na.rm = T)
tgc6$term <- "christian"
tgc7 <- summarySE(Cloud_MTurk, measurevar="vulnerable_eff", conf.interval = .90,na.rm = T)
tgc7$term <- "vulnerable"
tgc8 <- summarySE(Cloud_MTurk, measurevar="policy_eff", conf.interval = .90,na.rm = T)
tgc8$term <- "policy"
tgc9 <- summarySE(Cloud_MTurk, measurevar="econrisk_eff", conf.interval = .90,na.rm = T)
tgc9$term <- "economic risk"
tgc10 <- summarySE(Cloud_MTurk, measurevar="capacity_eff", conf.interval = .90,na.rm = T)
tgc10$term <- "capacity"
tgc11 <- summarySE(Cloud_MTurk, measurevar="muslim_eff", conf.interval = .90,na.rm = T)
tgc11$term <- "muslim"
tgc12 <- summarySE(Cloud_MTurk, measurevar="terrorist1_eff", conf.interval = .90,na.rm = T)
tgc12$term <- "terrorist 1"
tgc13 <- summarySE(Cloud_MTurk, measurevar="terrorist2_eff", conf.interval = .90,na.rm = T)
tgc13$term <- "terrorist 2"
tgc14 <- summarySE(Cloud_MTurk, measurevar="cultural_eff", conf.interval = .90,na.rm = T)
tgc14$term <- "cultural"
tgc15 <- summarySE(Cloud_MTurk, measurevar="aid_eff", conf.interval = .90,na.rm = T)
tgc15$term <- "aid"
tgc16 <- summarySE(Cloud_MTurk, measurevar="borders_eff", conf.interval = .90,na.rm = T)
tgc16$term <- "borders"
tgc17 <- summarySE(Cloud_MTurk, measurevar="disease_eff", conf.interval = .90,na.rm = T)
tgc17$term <- "disease"
tgc18 <- summarySE(Cloud_MTurk, measurevar="domore_eff", conf.interval = .90,na.rm = T)
tgc18$term <- "do more"
tgc19 <- summarySE(Cloud_MTurk, measurevar="christianreason_eff", conf.interval = .90,na.rm = T)
tgc19$term <- "christian reason"
tgc20 <- summarySE(Cloud_MTurk, measurevar="econ_eff", conf.interval = .90,na.rm = T)
tgc20$term <- "econ"


Cloud_MTurks <- c("tgc","tgc1","tgc2","tgc4","tgc5","tgc6","tgc7","tgc8","tgc9","tgc10","tgc11","tgc12","tgc13","tgc14","tgc15","tgc16","tgc17","tgc18","tgc19","tgc20")

for(tgcs in Cloud_MTurks)
  assign(tgcs, setNames(get(tgcs),  c("ID", "N", "mean","sd","se","ci","term")))

MASTER2 <- rbind(tgc,tgc1,tgc2,tgc4,tgc5,tgc6,tgc7,tgc8,tgc9,tgc10,tgc11,tgc12,tgc13,tgc14,tgc15,tgc16,tgc17,tgc18,tgc19,tgc20)

MASTER2$low = MASTER2$mean-MASTER2$ci
MASTER2$high = MASTER2$mean+MASTER2$ci

MASTER %>%
  left_join(MASTER2,by = "term") -> useme

plot <- ggplot(data=useme, aes(mean.y,mean.x)) +
  geom_point() +
  steph_theme + 
  labs(title="Mean Perceived Effectiveness & Support of Frames", 
       y="Very Opposed (1) - Very Supportive (7)", x = "Not very effective (1) - Very effective (4)", caption = "Winter 2020 CloudResearch + MTurk",subtitle = "Every Frame Tested:", fill = "Treatment Condition:")  +
  theme(legend.position = "none")

library(ggrepel)
plot + geom_text_repel(aes(label = useme$term),
                                size = 3.75) 


# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 1. Weighted Sample Demographics
# ---------------------------------------------------------------------------------------------------------------
questionr::wtd.table(x = df$gender_recode, weights = df$weight_bornagain,na.show = T)
weights::wpct(x = df$gender_recode, weight = df$weight_bornagain,na.rm = T)

questionr::wtd.table(x = df$educ_recode, weights = df$weight_bornagain,na.show = T)
weights::wpct(x = df$educ_recode, weight = df$weight_bornagain,na.rm = T)

questionr::wtd.table(x = df$faminc_new, weights = df$weight_bornagain,na.show = T)
100* weights::wpct(x = df$faminc_new, weight = df$weight_bornagain,na.rm = T)

questionr::wtd.table(x = df$ideo5_recode_num, weights = df$weight_bornagain,na.show = T)
100* weights::wpct(x = df$ideo5_recode_num, weight = df$weight_bornagain,na.rm = T)

df$age_recode <- car::recode(as.numeric(as.character(df$age)), "0:17 = '0:17'; 18:24 = '18:24'; 25:39 = '25:39'; 40:60 = '40:60'; 61:200 = '61+'; else = NA")
questionr::wtd.table(x = df$age_recode, weights = df$weight_bornagain,na.show = T)
100* weights::wpct(x = df$age_recode, weight = df$weight_bornagain,na.rm = T)

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 2. Balance Checks displaying multinomial logistic regression p-values.
# ---------------------------------------------------------------------------------------------------------------
regression <- multinom(messages ~ gender + age + faminc_new + educ + ideo5_recode_num, data = df, na.rm = T,weights = weight_bornagain,Hess = T) # Run regression
summary <- summary(regression)$coefficients/summary(regression)$standard.errors #store summary
p1 <- (1 - pnorm(abs(summary), 0, 1)) * 2 #two tail z test

df %>%
  mutate(messages = fct_relevel(df$Exp3_recode, "RV Message")) -> df_H2
regression <- multinom(messages ~ gender + age + faminc_new + educ + ideo5_recode_num, data = df_H2, na.rm = T,weights = weight_bornagain,Hess = T) # Run regression
summary <- summary(regression)$coefficients/summary(regression)$standard.errors #store summary
p2 <- (1 - pnorm(abs(summary), 0, 1)) * 2 #two tail z test

p1 <- as.data.frame(p1)
p2 <- as.data.frame(p2)
rownames(p1) <- c("RV+SC - Control", "RV - Control")
rownames(p2) <- c("Control - RV", "RV+SC - RV")
p2 <- p2[2,]
rbind(p1,p2)

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 3. Effect of the treatments on other groups.
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages, data = df, weights = weight_bornagain)
mod2 <- lm(FT_Evangelicals ~ messages, data = df, weights = weight_bornagain)
mod3 <- lm(FT_Republicans ~ messages, data = df, weights = weight_bornagain)
mod4 <- lm(FT_Democrats ~ messages, data = df, weights = weight_bornagain)
huxtable::huxreg('FT: Refugees' = mod1,
                 'FT: Evangelicals' = mod2, 
                 'FT: Republicans' = mod3, 
                 'FT: Democrats' = mod4, 
                 stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 4. Effect of the Religious Values + Source Cue message compared to the Religious Values message:
# ---------------------------------------------------------------------------------------------------------------
df %>%
  mutate(messages = fct_relevel(df$messages, "RV Message")) -> df_H2

mod1 <- lm(FT1_1_recode ~ messages, data = df_H2, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages, data = df_H2, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages, data = df_H2, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages, data = df_H2, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))


# ---------------------------------------------------------------------------------------------------------------
# Table 5. Treatment effects moderated by evangelical identity strength
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain)
mod2 <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain)
mod3 <- lm(BR1_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain)
mod4 <- lm(BR2_recode ~ messages * Evangelical_Index + Republican_Index, data = df, weights = weight_bornagain)
export_summs(mod1,mod2, mod3,mod4,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Table 6. Treatment effects moderated by Republican identity strength
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain)
mod2 <- lm(FuturePolicy_Scale ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain)
mod3 <- lm(BR1_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain)
mod4 <- lm(BR2_recode ~ messages * Republican_Index + Evangelical_Index, data = df, weights = weight_bornagain)
export_summs(mod1,mod2, mod3, mod4,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 7. Effect of the treatments moderated by voting for Trump.
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Figure 2. Effect of the treatments moderated by voting for Trump
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain,subset = messages != "Religious Values + Source Cue")
mod3c <- lm(FuturePolicy_Scale ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain,subset = messages != "Religious Values + Source Cue")
mod2c <- lm(BR1_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain,subset = messages != "Religious Values + Source Cue")
mod4c <- lm(BR2_recode ~ messages * E6_recode2, data = df_vote2, weights = weight_bornagain,subset = messages != "Religious Values + Source Cue")

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "E6_recode2", facet_labs = c("Religious Values + Source Cue * Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = No Trump Vote to 1 = Trump Vote", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=1)) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "E6_recode2", facet_labs = c("Religious Values + Source Cue * Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = No Trump Vote to 1 = Trump Vote", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=1)) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "E6_recode2", facet_labs = c("Religious Values + Source Cue * Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = No Trump Vote to 1 = Trump Vote", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=1)) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "E6_recode2", facet_labs = c("Religious Values + Source Cue * Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = No Trump Vote to 1 = Trump Vote", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=1)) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "E6_recode2", facet_labs = c("Religious Values * Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = No Trump Vote to 1 = Trump Vote", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=1)) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "E6_recode2", facet_labs = c("Religious Values * Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = No Trump Vote to 1 = Trump Vote", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=1)) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "E6_recode2", facet_labs = c("Religious Values * Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Did not Vote for Trumo to 1 = Voted for Trump", y = "Average Effect", caption = "Marginal effects of the Religious values message relative to the control group by the respondent's vote choice. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith evangelical identification. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=1)) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "E6_recode2", facet_labs = c("Religious Values * Unemployment"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = No Trump Vote to 1 = Trump Vote", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=1)) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4)


# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 8. Effect of the treatments moderated by Trump approval.
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * C3_recode, data = df_approve, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Figure 3. Effect of the treatments moderated by Trump approval
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * C3_recode, data = df_approve, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod2c <- lm(BR1_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod4c <- lm(BR2_recode ~ messages * C3_recode, data = df_approve, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod3c <- lm(FuturePolicy_Scale ~ messages * C3_recode, data = df_approve, weights = weight_bornagain,subset = messages != "RV + SC Message")
scaleFUN <- function(x) sprintf("%.3f", x)

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "C3_recode", facet_labs = c("Religious Values + Source Cue * Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Strongly Disapprove to 4 = Strongly Approve", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "C3_recode", facet_labs = c("Religious Values + Source Cue * Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Strongly Disapprove to 4 = Strongly Approve", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "C3_recode", facet_labs = c("Religious Values + Source Cue * Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Strongly Disapprove to 4 = Strongly Approve", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "C3_recode", facet_labs = c("Religious Values + Source Cue * Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Strongly Disapprove to 4 = Strongly Approve", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "C3_recode", facet_labs = c("Religious Values * Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Strongly Disapprove to 4 = Strongly Approve", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "C3_recode", facet_labs = c("Religious Values * Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Strongly Disapprove to 4 = Strongly Approve", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "C3_recode", facet_labs = c("Religious Values * Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Strongly Disapprove to 4 = Strongly Approve", y = "Average Effect", caption = "Marginal effects of the Religious values message relative to the control group by the respondent's support of Trump. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith support for Trump. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "C3_recode", facet_labs = c("Religious Values * Unemployment"),
                    point = T,ci = .90, plot = T) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Strongly Disapprove to 4 = Strongly Approve", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4)

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 9. Direct effects of the Religious Values and Religious Values + Source Cue messages on attitudes and 
# policy opinion (including speeders)
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 10. Effect of the treatments moderated by evangelical identity controlling for Republican strength (including speeders)
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * Evangelical_Index + Republican_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * Evangelical_Index + Republican_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * Evangelical_Index + Republican_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index + Republican_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Figure 4. Effect of the treatments moderated by evangelical identity with Republican controls (including speeders)
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod3c <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod2c <- lm(BR1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod4c <- lm(BR2_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
scaleFUN <- function(x) sprintf("%.3f", x)

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values + Source Cue * Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values + Source Cue * Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values + Source Cue * Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values + Source Cue * Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

scaleFUN <- function(x) sprintf("%.3f", x)

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values * Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))


plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values * Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))


plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values * Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "Marginal effects of the Religious values message relative to the control group by the respondent's evangelical identification. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith evangelical identification. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values * Unemployment"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4)

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 11. Effect of the treatments moderated by Republican identity with evangelical controls (including speeders)
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Figure 5. Effect of the treatments moderated by Republican identity with evangelical controls (including speeders)
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod3c <- lm(FuturePolicy_Scale ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod2c <- lm(BR1_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod4c <- lm(BR2_recode ~ messages * Republican_Index + Evangelical_Index, data = df_SPEEDERSINCLUDED, weights = weight_bornagain,subset = messages != "RV + SC Message")
scaleFUN <- function(x) sprintf("%.3f", x)

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values + Source Cue * Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values + Source Cue * Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values + Source Cue * Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values + Source Cue * Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

scaleFUN <- function(x) sprintf("%.3f", x)

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values * Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))


plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values * Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))


plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values * Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "Marginal effects of the Religious values message relative to the control group by the respondent's Republican identification. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith Republican identification. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values * Unemployment"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4)

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 12. Effects of the COVID in refugee camps experiment on all DVs
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ bcv_treat_recode, data = df, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ bcv_treat_recode, data = df, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ bcv_treat_recode, data = df, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ bcv_treat_recode, data = df, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 13. Effects of the treatment messages controlling for COVID in refugee camps experiment on all DVs
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages + bcv_treat_recode, data = df, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages + bcv_treat_recode, data = df, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages + bcv_treat_recode, data = df, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages + bcv_treat_recode, data = df, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))


# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 14. Effect of the treatments moderated by an alternate measure of evangelical identity
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * R3_recode, data = df, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * R3_recode, data = df, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * R3_recode, data = df, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * R3_recode, data = df, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Figure 6. Effect of the treatments moderated by an alternate measure of evangelical identity
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * R3_recode, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * R3_recode, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * R3_recode, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * R3_recode, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * R3_recode, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod3c <- lm(FuturePolicy_Scale ~ messages * R3_recode, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod2c <- lm(BR1_recode ~ messages * R3_recode, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod4c <- lm(BR2_recode ~ messages * R3_recode, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "R3_recode", facet_labs = c("Religious Values + Source Cue * Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Weak Evangelical to 4 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))


plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "R3_recode", facet_labs = c("Religious Values + Source Cue * Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Weak Evangelical to 4 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "R3_recode", facet_labs = c("Religious Values + Source Cue * Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Weak Evangelical to 4 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "R3_recode", facet_labs = c("Religious Values + Source Cue * Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Weak Evangelical to 4 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "R3_recode", facet_labs = c("Religious Values * Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Weak Evangelical to 4 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "R3_recode", facet_labs = c("Religious Values * Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Weak Evangelical to 4 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "R3_recode", facet_labs = c("Religious Values * Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Weak Evangelical to 4 = Strong Evangelical", y = "Average Effect", caption = "Marginal effects of the Religious values message relative to the control group by the respondent's evangelical identification. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith evangelical identification. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "R3_recode", facet_labs = c("Religious Values * Unemployment"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="1 = Weak Evangelical to 4 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4) 

# ---------------------------------------------------------------------------------------------------------------
# Appendix Figure 7. Religious message effects relative to the control group, by evangelical strength
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod3c <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod2c <- lm(BR1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod4c <- lm(BR2_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
scaleFUN <- function(x) sprintf("%.3f", x)

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values + Source Cue * Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values + Source Cue * Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values + Source Cue * Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values + Source Cue * Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

scaleFUN <- function(x) sprintf("%.3f", x)

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values * Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values * Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values * Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "Marginal effects of the Religious message relative to the control group by the respondent's Evangelical identification. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith Evangelical identification. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "Evangelical_Index", facet_labs = c("Religious Values * Unemployment"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Evangelical to 1 = Strong Evangelical", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4)

# ---------------------------------------------------------------------------------------------------------------
#Appendix Table 15. Effect of the treatments moderated by evangelical identity without Republican identity controls
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * Evangelical_Index, data = df, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index, data = df, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Figure 8. Religious message effects relative to the control group, by Republican strength without evangelical identity controls
# ---------------------------------------------------------------------------------------------------------------
mod1a <- lm(FT1_1_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod3a <- lm(FuturePolicy_Scale ~ messages * Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod2a <- lm(BR1_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod4a <- lm(BR2_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV Message")
mod1c <- lm(FT1_1_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod3c <- lm(FuturePolicy_Scale ~ messages * Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod2c <- lm(BR1_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
mod4c <- lm(BR2_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain,subset = messages != "RV + SC Message")
scaleFUN <- function(x) sprintf("%.3f", x)

plot1 <- interplot(m = mod1a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values + Source Cue * Feeling Thermometer"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2 <- interplot(m = mod3a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values + Source Cue * Resettlement"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3 <- interplot(m = mod2a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values + Source Cue * Schooling"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4 <- interplot(m = mod4a, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values + Source Cue * Unemployment"),
                   point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

scaleFUN <- function(x) sprintf("%.3f", x)

plot1b <- interplot(m = mod1c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values * Feeling Thermometer"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1), "cm"))

plot2b <- interplot(m = mod3c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values * Resettlement"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(0,1,-1,1.5), "cm"))

plot3b <- interplot(m = mod2c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values * Schooling"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "Marginal effects of the Religious message relative to the control group by the respondent's Republican identification. 90% confidence \nintervals (two-tailed) are illustrated. The marginal effects are from a set of analyses in which the treatment conditions were interacted \nwith Republican identification. OLS was used for all analyses.",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,0.5,1), "cm"))

plot4b <- interplot(m = mod4c, var1 = "messages", var2 = "Republican_Index", facet_labs = c("Religious Values * Unemployment"),
                    point = T,ci = .90) +
  steph_theme +
  theme(axis.text.x  = element_text(angle=90)) +
  geom_hline(yintercept = 0, linetype = "dashed") + 
  labs(title="", 
       x="0 = Weak Republican to 1 = Strong Republican", y = "Average Effect", caption = "",subtitle = "", fill = ":") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  scale_x_continuous(labels=scaleFUN) +
  theme(plot.margin=unit(c(-1,1,1,1.5), "cm"))

ggpubr::ggarrange(plot1,plot2,plot3,plot4,plot1b,plot2b,plot3b,plot4b, ncol = 2, nrow = 4)

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 16. Effect of the treatments moderated by Republican identity without evangelical identity controls
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * Republican_Index, data = df, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * Republican_Index, data = df, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 17. Direct effects of the religious Values and Religious Values + Source Cue messages on attitudes and policy opinion with demographic controls
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 18. Effect of the treatments moderated by evangelical identity with Republican strength and demographic controls
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * Evangelical_Index + Republican_Index +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * Evangelical_Index + Republican_Index  +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * Evangelical_Index + Republican_Index  +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * Evangelical_Index + Republican_Index  +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 19. Effect of the treatments moderated by Republican identity with evangelical strength and demographic controls
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * Republican_Index + Evangelical_Index +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod2 <- lm(BR1_recode ~ messages * Republican_Index + Evangelical_Index  +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod3 <- lm(BR2_recode ~ messages * Republican_Index + Evangelical_Index  +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
mod6 <- lm(FuturePolicy_Scale ~ messages * Republican_Index + Evangelical_Index  +  gender + age + faminc_new + educ, data = df, weights = weight_bornagain)
export_summs(mod1,mod6, mod2,mod3,
             error_format = "({std.error})",
             error_pos = c("below"),
             model.names = c('Feeling Thermometer','Resettlement','Schooling','Unemployment'),
             coefs = NULL, stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01))

# ---------------------------------------------------------------------------------------------------------------
# Appendix Table 20. Religious values message treatment effects in the CES among white non-Evangelicals
# ---------------------------------------------------------------------------------------------------------------
CCES %>%
  mutate(race_r = as_factor(race)) %>%
  mutate(pew_bornagain_r = as_factor(pew_bornagain)) %>%
  filter(pew_bornagain_r == "No",
         race_r == "White") -> df_nonevan_white

mod1 <- lm(FT1_1_recode ~ messages, data = df_nonevan_white, weights = teamweight)
mod2 <- lm(BR1_recode ~ messages, data = df_nonevan_white, weights = teamweight)
mod3 <- lm(BR2_recode ~ messages, data = df_nonevan_white, weights = teamweight)
mod6 <- lm(FuturePolicy_Scale ~ messages, data = df_nonevan_white, weights = teamweight)

export_summs(mod1,mod6,mod2,mod3,
             error_pos = c("below"),
             ci_level = 0.90,
             stars = c(`*` = 0.1, `**` = 0.05, `***` = 0.01),
             model.names = c('Feeling Therm.', 'Resettlement', 'Schooling','Unemployment'),
             coefs = NULL)


# ---------------------------------------------------------------------------------------------------------------
# Figure B1: Church Attendance in our Sample
# ---------------------------------------------------------------------------------------------------------------
ggplot(data = df, aes(x = pew_churatd_r), weight = weight_bornagain) +
  geom_histogram(binwidth = 0.16) +
  labs(title="Distibution of Church Attendance", 
       x="Low to High Attendance", caption = "Weighted YouGov Data",subtitle = "", fill = "Key:") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  steph_theme

# ---------------------------------------------------------------------------------------------------------------
# Figure B2: Frequency of Prayer in our Sample
# ---------------------------------------------------------------------------------------------------------------
ggplot(data = df, aes(x = pew_prayer_r), weight = weight_bornagain) +
  geom_histogram(binwidth = 0.16) +
  labs(title="Distibution of Frequency of Prayer", 
       x="Low to High Frequency", caption = "Weighted YouGov Data",subtitle = "", fill = "Key:") + 
  theme(legend.position = "none",axis.text.x = element_text(size=12, angle=35)) +
  steph_theme

# ---------------------------------------------------------------------------------------------------------------
# Table B1. Effect of the treatments moderated by church attendance
# ---------------------------------------------------------------------------------------------------------------
mod1 <- lm(FT1_1_recode ~ messages * pew_churatd_r, data = df, weights = weight_bornagain)
mod2 <- lm(FuturePolicy_Scale ~ messages * pew_churatd_r, data = df, weights = weight_bornagain)
mod3 <- lm(BR1_recode ~ messages * pew_churatd_r, data = df, weights = weight_bornagain)
mod4 <- lm(BR2_recode ~ messages * pew_churatd_r, data = df, weights = weight_bornagain)
export_summs(mod1,mod2,mod3,mod4,
             error_format = "({std.error})",
             error_pos = c("below"),
             ci_level = 0.90,
             model.names = c('Feeling Thermometer','Resettlement', 'Stay in School','Unemployment'),
             coefs = NULL)



